clear; close all; clc;

format long

% Note: the variables v1, v2, v3 are used to determine which regions have
% solid lines versus dashed lines (decreasing vs increasing parts of the
% schedule).

%% Benchmark case

cases = 1;
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 6: Interest rate spreads for Argentina (b = 20%)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=2;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>11;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<11)|(yy>20);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=2;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',2.00,'LineStyle','--'); hold on
ig=2;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
ig=2;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',2.00,'LineStyle','--');
ig=2;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',10)
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',10)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',10)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,'low growth','high growth');
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.5]; % [width height]
paper_position = [0 0 5.0 2.5]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rcorresp.pdf';
saveas(gcf,filename)

% Figure 7.b : Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Inherited debt b
bb1 = 0.20;
bb2 = 0.15;
ib1v  = find(bvec<bb1,1,'last');
ib2v  = find(bvec<bb2,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib2v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib2v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib2v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib2v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib2v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$b$=',num2str(round(100*bb2,0)),'\%'], ...
             ['$b$=',num2str(round(100*bb1,0)),'\%']);
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.0]; % [width height]
paper_position = [0 0 5.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rsch_b.pdf';
saveas(gcf,filename)

% Figure 10.a: Policy functions and equilibrium spreads for Argentina - Debt issuance (n), low growth
figure
ig=1;is=2;iz=floor((Nz+1)/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p1=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50); hold on
ig=1;is=1;iz=floor((Nz+1)/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p2=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 40])
ylim([10 40])
hleglines = [p1(1) p2(1)];
leg = legend(hleglines,'good sunspot','bad sunspot');
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_npol_gL.pdf';
saveas(gcf,filename)

% Figure 10.b: Policy functions and equilibrium spreads for Argentina - Debt issuance (n), high growth
figure
ig=2;is=2;iz=floor((Nz+1)/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p3=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50); hold on
ig=2;is=1;iz=floor((Nz+1)/2);xx=100*bvec;yy=100*nopt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p4=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 40])
ylim([10 40])
hleglines = [p3(1) p4(1)];
leg = legend(hleglines,'good sunspot','bad sunspot');
set(leg,'Interpreter','LaTex','Fontsize',8,'Location','SouthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_npol_gH.pdf';
saveas(gcf,filename)

% Figure 10.c: Policy functions and equilibrium spreads for Argentina - Spreads (ρ − r∗), low growth
figure
ig=1;is=2;iz=floor((Nz+1)/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p1=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50); hold on
ig=1;is=1;iz=floor((Nz+1)/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p2=plot(xx(cc),yy(cc),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9)
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^{*}$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 40])
ylim([0 40])
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rpol_gL.pdf';
saveas(gcf,filename)

% Figure 10.d: Policy functions and equilibrium spreads for Argentina - Spreads (ρ − r∗), high growth
figure
ig=2;is=2;iz=floor((Nz+1)/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p3=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50); hold on
ig=2;is=1;iz=floor((Nz+1)/2);xx=100*bvec;yy=Ropt(:,ig,is,iz,1);cc=(defpol(:,ig,is,iz,1)<0.5); p4=plot(xx(cc),yy(cc),'color',[0 0 0.8],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9)
xlabel('debt payments $b$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^{*}$ (\%)','Interpreter','LaTex','Fontsize',9)
xlim([0 40])
ylim([0 40])
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.0 2.0]; % [width height]
paper_position = [0 0 3.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_benchmark_Rpol_gH.pdf';
saveas(gcf,filename)

%% Varying the probability of the bad sunspot: psB

cases = [2 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 7.a : Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - The role of expectation pB
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>5;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<10)|(yy>30);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p5(1) p1(1)];
leg = legend(hleglines,...
             ['$p_B$=',num2str(round(param(12,1),2))], ...
             ['$p_B$=',num2str(round(param(12,2),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.0]; % [width height]
paper_position = [0 0 5.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_pBs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pBs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table 3: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Case 1: pB=',num2str(round(param(12,1),2)),'\n'])
fprintf(['Case 2: pB=',num2str(round(param(12,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Case 1               \t Case 2                \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying maturity (delta)

cases = [3 1 4];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);
ibv= zeros(ncases,1);
DD = 0.50;

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);    
    load(['OUTPUT_',num2str(cases(ic)),'/bvec.txt'])     
    ibv(ic)  = find(bvec<DD*param(2,ic),1,'last');             
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure 7.c : Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Maturity δ (b/δ = 50%)
figure
ig=1;is=1; xx = 100*nsch(:,ibv(1),ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>20;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<20)|(yy>35);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ibv(1),ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',3.00,'LineStyle',':'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ibv(1),ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',1.0,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ibv(1),ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',3.00,'LineStyle',':');
ig=1;is=1;p4=plot(100*nsch(v3:end,ibv(1),ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0.5 0],'LineWidth',1.0,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ibv(2),ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ibv(2),ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ibv(2),ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ibv(2),ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p8=plot(100*nsch(v3:end,ibv(2),ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ibv(3),ig,is,3); yy = 100*(Rsch(:,ig,is,3)-(param(1,3)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p9=plot(100*nsch(1:v1,ibv(3),ig,is,3),100*(Rsch(1:v1,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p10=plot(100*nsch(v1:v2,ibv(3),ig,is,3),100*(Rsch(v1:v2,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p11=plot(100*nsch(v2:v3,ibv(3),ig,is,3),100*(Rsch(v2:v3,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p12=plot(100*nsch(v3:end,ibv(3),ig,is,3),100*(Rsch(v3:end,ig,is,3)-(param(1,3)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1) p9(1)];
leg = legend(hleglines,...
             ['$\delta$=',num2str(round(param(2,1),2)),' (' ,num2str(round(1/param(2,1),1)),' years)'], ...
             ['$\delta$=',num2str(round(param(2,2),2)),' (' ,num2str(round(1/param(2,2),1)),' years)'], ...
             ['$\delta$=',num2str(round(param(2,3),2)),' (' ,num2str(round(1/param(2,3),1)),' years)']); 
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.0]; % [width height]
paper_position = [0 0 5.0 2.0]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_deltas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_deltas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (2): delta=',num2str(round(param(2,1),2)),'\n'])
fprintf(['Column (1): delta=',num2str(round(param(2,2),2)),'  (benchmark) \n'])
fprintf(['Column (3): delta=',num2str(round(param(2,3),2)),'\n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (2)           \t Column (1)           \t Column (3)            \n')
fprintf('-------------------                \t -------------------- \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying persistence of low-growth state (pL)

cases = [5 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);              
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7a: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Probability of remaining in low growth (pL)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$p_{L}$=',num2str(round(param(6,2),2))], ...
             ['$p_{L}$=',num2str(round(param(6,1),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_pgLs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pgLs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (4): pgLs=',num2str(round(param(6,1),2)),'\n'])
fprintf(['Column (1): pgLs=',num2str(round(param(6,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (4)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the recovery rate (kappa)

cases = [6 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);               
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7b: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Recovery rate (κ)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$\kappa$=',num2str(round(param(5,2),2))], ...
             ['$\kappa$=',num2str(round(param(5,1),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_kappas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_kappas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
         
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (5): kappa=',num2str(round(param(5,1),2)),'\n'])
fprintf(['Column (1): kappa=',num2str(round(param(5,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (5)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the dispersion of idiosyncratic shocks (sigma)

cases = [7 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                 
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7c: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - St. deviation of transitory shock (σ)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',3.00,'LineStyle',':'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.0,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',3.00,'LineStyle',':');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.0,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$\sigma$=',num2str(round(100*param(11,2),1)),'\%'], ...
             ['$\sigma$=',num2str(round(100*param(11,1),1)),'\%']); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_sigmas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_sigmas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (6): sigma=',num2str(round(param(11,1),3)),'\n'])
fprintf(['Column (1): sigma=',num2str(round(param(11,2),3)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (6)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the default cost in the high-growth state (ybarH)

cases = [8 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7d: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Cost of default in high growth (ϕ(gH))
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$\phi (g_{H})$=',num2str(round(param(4,2),2))], ...
             ['$\phi (g_{H})$=',num2str(round(param(4,1),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_ybarHs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_ybarHs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (7): ybarH=',num2str(round(param(4,1),2)),'\n'])
fprintf(['Column (1): ybarH=',num2str(round(param(4,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (7)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the persistence of the high-growth stat (pH)

cases = [9 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                 
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7e: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Probability of remaining in high growth (pH)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$p_{H}$=',num2str(round(param(7,2),2))], ...
             ['$p_{H}$=',num2str(round(param(7,1),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_pgHs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_pgHs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (8): pH=',num2str(round(param(7,1),2)),'\n'])
fprintf(['Column (1): pH=',num2str(round(param(7,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (8)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the discount factor (beta)

cases = [10 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7f: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Discount factor (β)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$\beta$=',num2str(round(param(8,2),2))], ...
             ['$\beta$=',num2str(round(param(8,1),2))]); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_betas_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_betas.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (9): beta=',num2str(round(param(8,1),2)),'\n'])
fprintf(['Column (1): beta=',num2str(round(param(8,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (9)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Varying the difference between high and low growth rates (gH-gL)

cases = [11 1];
ncases = length(cases);
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])

% Parameters
Tburn   = 1000;
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
Rsch = zeros(Nb,Ng,Ns,ncases);
nsch = zeros(Nb,Nb,Ng,Ns,ncases);
defpol = zeros(Nb,Ng,Ns,Nz,ncases);
probdef     = zeros(Nb,Ng,Ns,ncases);
pol_probdef = zeros(Nb,Ng,Ns,Nz,ncases);
bopt   = zeros(Nb,Ng,Ns,Nz,ncases);
bbopt  = zeros(Nb,Ng,Ns,Nz,ncases);
nopt   = zeros(Nb,Ng,Ns,Nz,ncases);
Ropt   = zeros(Nb,Ng,Ns,Nz,ncases);
% Moments
avg_rho = zeros(ncases,1); avg_rho_gL = zeros(ncases,1); avg_rho_gH = zeros(ncases,1);
avg_qb = zeros(ncases,1);avg_qb_gL = zeros(ncases,1); avg_qb_gH = zeros(ncases,1);
avg_cn = zeros(ncases,1); avg_cn_gL = zeros(ncases,1); avg_cn_gH = zeros(ncases,1);
avg_n = zeros(ncases,1); avg_n_gL = zeros(ncases,1); avg_n_gH = zeros(ncases,1);
avg_b = zeros(ncases,1); avg_b_gL = zeros(ncases,1); avg_b_gH = zeros(ncases,1);
avg_tb = zeros(ncases,1); avg_tb_gL = zeros(ncases,1); avg_tb_gH = zeros(ncases,1);
def_rate = zeros(ncases,1); def_rate_gL = zeros(ncases,1); def_rate_gH = zeros(ncases,1); 
std_c_y = zeros(ncases,1);
corr_rho_y = zeros(ncases,1);
std_rho = zeros(ncases,1);

% prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
param  = zeros(Nparam,ncases);
for ic = 1:ncases
    % Parameters
    load(['OUTPUT_',num2str(cases(ic)),'/prms.txt'])    
    for j = 1:Nparam
        param(j,ic) = prms(j);
    end
    dlta = param(2,ic);                 
    % Schedule
    load(['OUTPUT_',num2str(cases(ic)),'/Rsched.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/Ednp.txt'])        
    count = 0;
    for ig = 1:Ng
        for is = 1:Ns
            for ib = 1:Nb
                count = count+1;
                Rsch(ib,ig,is,ic) = Rsched(count);
                probdef(ib,ig,is,ic) = Ednp(count);
            end
        end
    end   
    for ig = 1:Ng
    for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsch(:,ib,ig,is,ic) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsch(:,ig,is,ic));
    end
    end
    end
    Rsch(:,:,:,ic) = dlta*Rsch(:,:,:,ic)-dlta; 
    % Policy functions
    load(['OUTPUT_',num2str(cases(ic)),'/dpol.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/bopt_loc.txt'])
    count = 0;
    for ig = 1:Ng
    for is = 1:Ns
    for iz = 1:Nz
    for ib = 1:Nb
        count = count+1;
        defpol(ib,ig,is,iz,ic) = dpol(count);
        bopt(ib,ig,is,iz,ic)   = bopt_loc(count);
        bbopt(ib,ig,is,iz,ic)  = bvec(bopt(ib,ig,is,iz,ic));
        nopt(ib,ig,is,iz,ic)   = nsch(bopt(ib,ig,is,iz,ic),ib,ig,is,ic);
        Ropt(ib,ig,is,iz,ic)   = 100*(Rsch(bopt(ib,ig,is,iz,ic),ig,is,ic)-(param(1,ic)-1));
        pol_probdef(ib,ig,is,iz,ic) = 100*probdef(bopt(ib,ig,is,iz,ic),ig,is,ic);
    end
    end
    end 
    end             
    % Simulation moments
    load(['OUTPUT_',num2str(cases(ic)),'/bsim.txt']) 
    load(['OUTPUT_',num2str(cases(ic)),'/cnsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/idsim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/rhosim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/nsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/gsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/zsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/issim.txt'])
    load(['OUTPUT_',num2str(cases(ic)),'/ysim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qbsim.txt'])    
    load(['OUTPUT_',num2str(cases(ic)),'/qsim.txt'])            
    % Adjusting grids and defining default event
    idefsim = idsim*0;
    idefsim(2:end) = ((idsim(2:end)-idsim(1:end-1))>0.95);
    bpsim   = bsim(2:end);
    bsim    = bsim(1:end-1);    
    cnsim   = cnsim(1:end-1);
    idsim   = idsim(2:end);
    idefsim = idefsim(2:end);
    % Excluding first Tburn periods
    bsim    = bsim(Tburn+1:end);
    qbsim   = qbsim(Tburn+1:end);
    qsim    = qsim(Tburn+1:end);    
    cnsim   = cnsim(Tburn+1:end);
    idsim   = idsim(Tburn+1:end);
    idefsim = idefsim(Tburn+1:end);
    rhosim  = rhosim(Tburn+1:end);
    rhosim  = rhosim-(param(1,ic)-1);
    nsim    = nsim(Tburn+1:end);
    dysim   = log(gsim(Tburn+1:end))+zsim(Tburn+1:end)-zsim(Tburn:end-1);    
    gsim    = gsim(Tburn+1:end);
    zsim    = zsim(Tburn+1:end);
    issim   = issim(Tburn+1:end);
    ysim    = ysim(Tburn+1:end);
    % Computing levels and normalizing by current GDP
    bsim    = bsim./ysim;
    qbsim   = qbsim./ysim;    
    cnsim   = cnsim./ysim;
    nsim    = nsim./ysim;    
    csim    = 1+nsim-bsim;
    tbsim   = 1-csim;
    % Moments excluding default periods
    avg_rho(ic,1) = 100*mean(rhosim(idsim==0));
    avg_rho_gL(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==min(gsim)))); avg_rho_gH(ic,1) = 100*mean(rhosim((idsim==0)&(gsim==max(gsim))));
    avg_qb(ic,1) = 100*mean(qbsim(idsim==0));
    avg_qb_gL(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==min(gsim)))); avg_qb_gH(ic,1) = 100*mean(qbsim((idsim==0)&(gsim==max(gsim))));
    avg_cn(ic,1) = 100*mean(cnsim(idsim==0));
    avg_cn_gL(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==min(gsim)))); avg_cn_gH(ic,1) = 100*mean(cnsim((idsim==0)&(gsim==max(gsim))));
    avg_n(ic,1) = 100*mean(nsim(idsim==0));
    avg_n_gL(ic,1) = 100*mean(nsim((idsim==0)&(gsim==min(gsim)))); avg_n_gH(ic,1) = 100*mean(nsim((idsim==0)&(gsim==max(gsim))));
    avg_b(ic,1) = 100*mean(bsim(idsim==0));
    avg_b_gL(ic,1) = 100*mean(bsim((idsim==0)&(gsim==min(gsim)))); avg_b_gH(ic,1) = 100*mean(bsim((idsim==0)&(gsim==max(gsim))));
    avg_tb(ic,1) = 100*mean(tbsim(idsim==0));
    avg_tb_gL(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==min(gsim)))); avg_tb_gH(ic,1) = 100*mean(tbsim((idsim==0)&(gsim==max(gsim))));
    def_rate(ic,1) = 100*(mean(idefsim)); def_rate_gL(ic,1) = 100*(mean(idefsim(gsim==min(gsim)))); def_rate_gH(ic,1) = 100*(mean(idefsim(gsim==max(gsim))));   
    corr_rho_y2 = corrcoef(rhosim(idsim==0),ysim(idsim==0));
    corr_rho_y(ic,1) = corr_rho_y2(1,2);
    std_rho(ic,1) = 100*std(rhosim(idsim==0));
    std_c_y(ic,1) = std(csim(idsim==0).*ysim(idsim==0))/std(ysim(idsim==0));
end

% Figure F.7g: Interest rate spreads for Argentina in the low growth state
% (b = 20%): Comparative statics - Difference between growth rates (gH − gL)
bb = 0.20;
ib1v  = find(bvec<bb,1,'last');
figure
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,2); yy = 100*(Rsch(:,ig,is,2)-(param(1,2)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p1=plot(100*nsch(1:v1,ib1v,ig,is,2),100*(Rsch(1:v1,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-'); hold on
ig=1;is=1;p2=plot(100*nsch(v1:v2,ib1v,ig,is,2),100*(Rsch(v1:v2,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p3=plot(100*nsch(v2:v3,ib1v,ig,is,2),100*(Rsch(v2:v3,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',2.00,'LineStyle','-');
ig=1;is=1;p4=plot(100*nsch(v3:end,ib1v,ig,is,2),100*(Rsch(v3:end,ig,is,2)-(param(1,2)-1)),'color',[0.8 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1; xx = 100*nsch(:,ib1v,ig,is,1); yy = 100*(Rsch(:,ig,is,1)-(param(1,1)-1));
ind1 = yy>15;
xx_aux = xx+(ind1*(-1.0e10));
[aa v1] = max(xx_aux);
ind2 = (yy<15)|(yy>38);
xx_aux = xx+(ind2*(1.0e10));
[aa v2] = min(xx_aux);
xx_aux = xx+(ind2*(-1.0e10));
[aa v3] = max(xx_aux);
ig=1;is=1;p5=plot(100*nsch(1:v1,ib1v,ig,is,1),100*(Rsch(1:v1,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--'); hold on
ig=1;is=1;p6=plot(100*nsch(v1:v2,ib1v,ig,is,1),100*(Rsch(v1:v2,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
ig=1;is=1;p7=plot(100*nsch(v2:v3,ib1v,ig,is,1),100*(Rsch(v2:v3,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',2.00,'LineStyle','--');
ig=1;is=1;p8=plot(100*nsch(v3:end,ib1v,ig,is,1),100*(Rsch(v3:end,ig,is,1)-(param(1,1)-1)),'color',[0 0 0],'LineWidth',1.50,'LineStyle',':');
set(gca,'XGrid','off','YGrid','off','Fontsize',9) 
xlabel('debt issuance $n$ (\%)','Interpreter','LaTex','Fontsize',9)
ylabel('spread $\rho-r^*$ ($\%$)','Interpreter','LaTex','Fontsize',9)
xlim([0 45])
ylim([0 40])
hleglines = [p1(1) p5(1)];
leg = legend(hleglines,...
             ['$\Delta^{g}$=',num2str(round(100*(param(10,2)-param(9,2)),0)),'\%'], ...
             ['$\Delta^{g}$=',num2str(round(100*(param(10,1)-param(9,1)),0)),'\%']); 
set(leg,'Interpreter','LaTex','Fontsize',7,'Location','SouthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [3.2 1.8]; % [width height]
paper_position = [0 0 3.2 1.8]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_gHgLs_Rsch.pdf';
saveas(gcf,filename)

filename = 'Results/Table_gHgLs.txt';

if (exist(filename,'file')==2)
   delete(filename);
end

diary(filename)
        
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Table F.2: Simulation moments: Argentina \n')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf(['Column (10): gH-gL=',num2str(round(param(10,1)-param(9,1),2)),'\n'])
fprintf(['Column (1):  gH-gL=',num2str(round(param(10,2)-param(9,2),2)),'  (benchmark) \n'])
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('Statistic                          \t Column (10)           \t Column (1)            \n')
fprintf('-------------------                \t -------------------- \t --------------------  \n')
fprintf('First moments                      \t                      \t                       \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Low-growth state                   \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gL(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gL(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gL(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gL(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gL(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gL(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gL(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('High-growth state                  \t                      \t                      \n')
fprintf('avg(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',avg_rho_gH(:,1)')
fprintf('avg(qb/dlta)(%%GDP)                \t %5.0f                \t %5.0f                \n',avg_qb_gH(:,1)'./param(2,:))
fprintf('avg(f)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_cn_gH(:,1)')
fprintf('avg(n)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_n_gH(:,1)')
fprintf('avg(b)(%%GDP)                      \t %5.0f                \t %5.0f                \n',avg_b_gH(:,1)')
fprintf('avg(tb) (%%GDP)                    \t %5.1f                \t %5.1f                \n',avg_tb_gH(:,1)')
fprintf('def rate (%%)                      \t %5.1f                \t %5.1f                \n',def_rate_gH(:,1)')
fprintf('                                   \t                      \t                      \n')
fprintf('Second moments                     \t                      \t                      \n')
fprintf('corr_rho_y                         \t %5.2f                \t %5.2f                \n',corr_rho_y(:,1)')
fprintf('std(rho-r*) (%%)                   \t %5.1f                \t %5.1f                \n',std_rho(:,1)')
fprintf('std(c)/std(y)                      \t %5.1f                \t %5.1f                \n',std_c_y(:,1)')
fprintf('----------------------------------------------------------------------------------------------------------- \n')
fprintf('\n \n \n')

diary off

%% Event case study: Argentina's 2001 default

cases = 1;
dirname = ['OUTPUT_',num2str(cases(1))];
load(['OUTPUT_',num2str(cases(1)),'/prms.txt'])
load(['OUTPUT_',num2str(cases(1)),'/bvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/gvec.txt'])
load(['OUTPUT_',num2str(cases(1)),'/zvec.txt'])
Ng   = length(gvec);
Ns   = 2;
Nb   = length(bvec);
Nz   = length(zvec);
Nparam = length(prms);
bopt   = zeros(Nb,Ng,Ns,Nz);
r_pol  = zeros(Nb,Ng,Ns,Nz);
% Policy functions
load(['OUTPUT_',num2str(cases(1)),'/bopt_loc.txt'])
load(['OUTPUT_',num2str(cases(1)),'/rpol.txt'])
count = 0;
for ig = 1:Ng
    for is = 1:Ns
        for iz = 1:Nz
            for ib = 1:Nb
                count = count+1;
                bopt(ib,ig,is,iz)   = bopt_loc(count);
                r_pol(ib,ig,is,iz)  = rpol(count);
            end
        end
    end
end
rs = prms(1)-1;
years    = (1997:1:2001);
T = length(years);
dy_data      =    [7.8;3.8;-3.4;-0.8;-4.5]/100;
r_data       =    [4.61;7.07;5.33;7.73;33.72];
ig_data      =    [2,2,1,1,1];
g_data       =    [0.04;0.04;-0.04;-0.04;-0.04];
iz_sim       =    0*dy_data;
iLz_sim      =    0*dy_data;
z_sim        =    0*dy_data;
Lz_sim       =    0*dy_data;
dy_sim       =    0*dy_data;
iz_sim(1)    = 12;
z_sim(1)     = zvec(iz_sim(1));
[a,bb]       = min(abs(dy_data(1)-(g_data(1)+z_sim(1)-zvec)));
iLz_sim(1)   = bb;
Lz_sim(1)    = zvec(iLz_sim(1));
dy_sim(1)    = g_data(1)+z_sim(1)-Lz_sim(1);
for t=2:T
    iLz_sim(t) = iz_sim(t-1);
    Lz_sim(t)  = z_sim(t-1);
    [a,bb]     = min(abs(dy_data(t)-(g_data(t)-Lz_sim(t)+zvec)));
    iz_sim(t)  = bb;
    z_sim(t)   = zvec(iz_sim(t));
    dy_sim(t)  = g_data(t)+z_sim(t)-Lz_sim(t);
end
% Sample path [97 98 99 00 01];
ib0 = 340;
ib1_sim = [ib0 0 0 0 0];
ib2_sim = [ib0 0 0 0 0];
is1_sim  = [1 1 1 1 1];
is2_sim  = [1 1 1 1 2];
r1_sim   = ib1_sim*0;
r2_sim   = ib2_sim*0;
for t=1:length(ib1_sim)
    if t<length(ib1_sim)
        ib1_sim(t+1)=bopt(ib1_sim(t),ig_data(t),is1_sim(t),iz_sim(t));
        ib2_sim(t+1)=bopt(ib2_sim(t),ig_data(t),is2_sim(t),iz_sim(t));
    end
     r1_sim(t)  = 100*(r_pol(ib1_sim(t),ig_data(t),is1_sim(t),iz_sim(t))-rs);    
     r2_sim(t)  = 100*(r_pol(ib2_sim(t),ig_data(t),is2_sim(t),iz_sim(t))-rs);    
end

% Figure 12: Argentina: Spreads from 1997 to 2001
figure
p1 = plot(years,r1_sim,'Color',[0.8 0 0],'LineStyle','-','LineWidth',1); hold on
p2 = plot(years,r2_sim,'Color',[0 0 0.8],'LineStyle','--','LineWidth',1);
p3 = plot(years,r_data,'Color',[0 0 0],'LineStyle',':','LineWidth',1.5);
xlim([1997,2001])
ylim([0,40])
xticks([1997,1998,1999,2000,2001])
ylabel('spreads (\%)','Interpreter','LaTex','Fontsize',10)
set(gca,'XGrid','off','YGrid','off','Fontsize',10)
hleglines = [p1(1) p2(1) p3(1)];
leg = legend(hleglines,...
             ['sequence with bad sunspot in 2001'], ...
             ['sequence with good sunspot in 2001'], ...
             ['data']);
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthWest')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.5]; % [width height]
paper_position = [0 0 5.0 2.5]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_sim_spreads.pdf';
saveas(gcf,filename)

% Figure E.1: Argentina: GDP growth from 1997 to 2001
figure
p1 = plot(years,100*dy_sim,'Color',[0.8 0 0],'LineStyle','-','LineWidth',1); hold on
p2 = plot(years,100*dy_data,'Color',[0 0 0],'LineStyle','--','LineWidth',1);
xlim([1997,2001])
ylim([-8,10])
xticks([1997,1998,1999,2000,2001])
ylabel('GDP growth (\%)','Interpreter','LaTex','Fontsize',10)
hleglines = [p1(1) p2(1)];
leg = legend(hleglines,...
             ['simulation'], ...
             ['data']);
set(leg,'Interpreter','LaTex','Fontsize',9,'Location','NorthEast')
legend boxoff
% specification
paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
paper_size = [5.0 2.5]; % [width height]
paper_position = [0 0 5.0 2.5]; % [left bottom width height].
paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
fig = gcf;
fig.PaperUnits = paper_unit;
fig.PaperSize = paper_size;
fig.PaperPosition = paper_position;
fig.PaperOrientation = paper_orientation;
filename = 'Results/fig_sim_gdp.pdf';
saveas(gcf,filename)